if (!require(pacman)) install.packages("pacman")
## Loading required package: pacman
pacman::p_load(tidyverse,
               sf,
               mapview,
               here)


# read/export vector data -------------------------------------------------

# read a shapefile (e.g., ESRI Shapefile format)
# `quiet = TRUE` just for cleaner output
(sf_nc_county <- st_read(dsn = here("data/nc.shp"),
                         quiet = TRUE))
## Simple feature collection with 100 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
## Geodetic CRS:  WGS 84
## First 10 features:
##         county                       geometry
## 1         ashe MULTIPOLYGON (((-81.47258 3...
## 2    alleghany MULTIPOLYGON (((-81.23971 3...
## 3        surry MULTIPOLYGON (((-80.45614 3...
## 4    currituck MULTIPOLYGON (((-76.00863 3...
## 5  northampton MULTIPOLYGON (((-77.21736 3...
## 6     hertford MULTIPOLYGON (((-76.74474 3...
## 7       camden MULTIPOLYGON (((-76.00863 3...
## 8        gates MULTIPOLYGON (((-76.56218 3...
## 9       warren MULTIPOLYGON (((-78.30849 3...
## 10      stokes MULTIPOLYGON (((-80.02545 3...
# export as shp file 
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.shp"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `ESRI Shapefile'
## Writing layer `sf_nc_county' to data source 
##   `/Users/kayahamilton/Documents/GitHub/cw_rgis/data/sf_nc_county.shp' using driver `ESRI Shapefile'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
# export as geopackage 
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.gpkg"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `GPKG'
## Writing layer `sf_nc_county' to data source 
##   `/Users/kayahamilton/Documents/GitHub/cw_rgis/data/sf_nc_county.gpkg' using driver `GPKG'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
# export as rds 
saveRDS(sf_nc_county,
        file = here("data/sf_nc_county.rds"))

# read rds
sf_nc_county <- readRDS(file = here("data/sf_nc_county.rds"))

# point data
sf_site <- readRDS(file = here("data/sf_finsync_nc.rds"))
sf_site
## Simple feature collection with 122 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.02452 ymin: 34.53497 xmax: -76.74611 ymax: 36.54083
## Geodetic CRS:  WGS 84
## # A tibble: 122 × 2
##    site_id                          geometry
##  * <chr>                         <POINT [°]>
##  1 finsync_nrs_nc-10013 (-81.51025 36.11188)
##  2 finsync_nrs_nc-10014 (-80.35989 35.87616)
##  3 finsync_nrs_nc-10020 (-81.74472 35.64379)
##  4 finsync_nrs_nc-10023  (-82.77898 35.6822)
##  5 finsync_nrs_nc-10024 (-77.75384 35.38553)
##  6 finsync_nrs_nc-10027 (-83.69678 35.02467)
##  7 finsync_nrs_nc-10029 (-80.65668 35.16119)
##  8 finsync_nrs_nc-10034 (-82.04497 36.09917)
##  9 finsync_nrs_nc-10041 (-80.46558 35.04365)
## 10 finsync_nrs_nc-10049 (-77.96322 34.58249)
## # ℹ 112 more rows
mapview(sf_site,
        col.regions = "black", # point's fill color
        legend = FALSE) # disable legend
## take the first 10 sites
sf_sites10 <- sf_site %>% 
  slice(1:10)

mapview(sf_sites10,
        color = "black", # line's color
        legend = FALSE) # disable legend
#line data 
(sf_str <- readRDS(here("data/sf_stream_gi.rds")))
## Simple feature collection with 5012 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -80.04671 ymin: 35.9001 xmax: -79.53284 ymax: 36.25706
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-79.90748 36.18... fid000001
## 2  LINESTRING (-79.89768 36.20... fid000002
## 3  LINESTRING (-79.94756 36.15... fid000003
## 4  LINESTRING (-79.94152 36.25... fid000004
## 5  LINESTRING (-79.94206 36.20... fid000005
## 6  LINESTRING (-79.90003 36.18... fid000006
## 7  LINESTRING (-79.94737 36.20... fid000007
## 8  LINESTRING (-79.87595 36.18... fid000008
## 9  LINESTRING (-79.88022 36.25... fid000009
## 10 LINESTRING (-79.94859 36.25... fid000010
mapview(sf_str,
        color = "steelblue", # line's color
        legend = FALSE) # disable legend
## take the first 10 sites
sf_str10 <- sf_str %>% 
  slice(1:10)

mapview(sf_str10,
        color = "steelblue", # line's color
        legend = FALSE) # disable legend
#polygon
mapview(sf_nc_county,
        col.regions = "tomato",
        legend = FALSE)
 # pick guilford county 
sf_nc_gi <- sf_nc_county %>% 
  filter(county == "guilford")

mapview(sf_nc_gi,
        col.regions = "salmon",
        legend = FALSE)
# use ggplot to visualize a map 
## not a great map
ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_str) +
  geom_sf(data = sf_site)

## a little better
ggplot() +
  geom_sf(data = sf_nc_gi) +
  geom_sf(data = sf_str)

# exercise ----------------------------------------------------------------
#1
sf_str_as <- readRDS(file = here("data/sf_stream_as.rds"))
sf_str_as
## Simple feature collection with 6361 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -81.73891 ymin: 36.24445 xmax: -81.24562 ymax: 36.588
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-81.32158 36.47... fid000001
## 2  LINESTRING (-81.31133 36.47... fid000002
## 3  LINESTRING (-81.38594 36.45... fid000003
## 4  LINESTRING (-81.39307 36.46... fid000004
## 5  LINESTRING (-81.3454 36.533... fid000005
## 6  LINESTRING (-81.34719 36.53... fid000006
## 7  LINESTRING (-81.3328 36.518... fid000007
## 8  LINESTRING (-81.34964 36.50... fid000008
## 9  LINESTRING (-81.35949 36.53... fid000009
## 10 LINESTRING (-81.35094 36.52... fid000010
#2
print(sf_str_as)
## Simple feature collection with 6361 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -81.73891 ymin: 36.24445 xmax: -81.24562 ymax: 36.588
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-81.32158 36.47... fid000001
## 2  LINESTRING (-81.31133 36.47... fid000002
## 3  LINESTRING (-81.38594 36.45... fid000003
## 4  LINESTRING (-81.39307 36.46... fid000004
## 5  LINESTRING (-81.3454 36.533... fid000005
## 6  LINESTRING (-81.34719 36.53... fid000006
## 7  LINESTRING (-81.3328 36.518... fid000007
## 8  LINESTRING (-81.34964 36.50... fid000008
## 9  LINESTRING (-81.35949 36.53... fid000009
## 10 LINESTRING (-81.35094 36.52... fid000010
print(sf_nc_county)
## Simple feature collection with 100 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
## Geodetic CRS:  WGS 84
## First 10 features:
##         county                       geometry
## 1         ashe MULTIPOLYGON (((-81.47258 3...
## 2    alleghany MULTIPOLYGON (((-81.23971 3...
## 3        surry MULTIPOLYGON (((-80.45614 3...
## 4    currituck MULTIPOLYGON (((-76.00863 3...
## 5  northampton MULTIPOLYGON (((-77.21736 3...
## 6     hertford MULTIPOLYGON (((-76.74474 3...
## 7       camden MULTIPOLYGON (((-76.00863 3...
## 8        gates MULTIPOLYGON (((-76.56218 3...
## 9       warren MULTIPOLYGON (((-78.30849 3...
## 10      stokes MULTIPOLYGON (((-80.02545 3...
#sf_nc_as Geodetic CRS: WGS 84
#sf_nc_county  Geodetic CRS: WGS 84

#4 
ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_str_as) 

#5
sf_nc_as <- sf_nc_county %>% 
  filter(county == "ashe")

ggplot() +
  geom_sf(data = sf_nc_as) +
  geom_sf(data = sf_str_as)